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Abstract — Previously  we  showed  how  delay  communication 
between  globally  coupled  self-propelled  agents  causes  new 
spatio-temporal  patterns  to  arise  when  the  delay  coupling  is 
fixed  among  all  agents  [1].  In  this  paper,  we  show  how  discrete, 
randomly  distributed  delays  affect  the  dynamical  patterns.  In 
particular,  we  investigate  how  the  standard  deviation  of  the  time 
delay  distribution  affects  the  stability  of  the  different  patterns 
as  well  as  the  switching  probability  between  coherent  states. 

I.  INTRODUCTION 

Numerous  recent  investigations  have  been  devoted  to  the 
study  of  interacting  multi-agent  or  swarming  systems  in 
various  natural  and  engineering  fields  of  study.  Investiga¬ 
tions  of  interacting  systems  have  revealed  the  emergence 
of  highly  complex  dynamic  behaviors  in  space  and  time 
which  arise  even  though  the  dynamics  of  a  single  agent  is 
quite  simple  [2].  In  particular,  these  multi-agent  swarms  can 
self-organize  in  complicated  spatio-temporal  patterns  that 
depend  on  the  details  of  the  inter-agent  interactions.  These 
investigations  have  been  motivated  by  and  had  an  impact  on 
many  diverse  biological  systems  such  as  bacterial  colonies, 
schooling  fish,  flocking  birds,  swarming  locusts,  ants,  and 
pedestrians  [3],  [4],  [5],  [6],  [7].  In  this  paper,  we  are 
interested  in  the  application  that  biological  analogies  have  on 
the  design  of  systems  of  autonomous,  inter-communicating 
robotic  systems  [8],  [9],  [10],  [11]  and  mobile  sensor  net¬ 
works  [12]. 

There  is  great  interest  to  design  agent-interaction  protocols 
to  carry  out  robotic  motion  planning,  consensus  and  cooper¬ 
ative  control,  and  spatio-temporal  formation.  One  methodol¬ 
ogy  is  to  combine  inter-agent  potentials  with  external  ones  in 
order  to  achieve  multi-agent  cooperative  motion  in  a  manner 
that  is  not  too  sensitive  with  respect  the  number  of  agents. 
Some  important  applications  making  use  of  scalable  numbers 
of  agents  are:  obstacle  avoidance  [10],  boundary  tracking 
[13],  [14],  environmental  sensing  [12],  [15],  decentralized 
target  tracking  [16],  environmental  consensus  estimation 
[12],  [17]  and  task  allocation  [18], 

Authors  have  employed  very  diverse  approaches  in  the 
study  of  multi-agent  systems.  Some  authors  have  described 
the  swarms  at  the  individual  level,  writing  their  models  in 
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terms  of  ordinary  differential  equations  (ODEs)  or  delay 
differential  equations  (DDEs)  to  describe  their  trajectories 
[19],  [20],  [9].  The  addition  of  noise  on  the  swarm’s  dynam¬ 
ics  introduces  even  richer  behavior,  such  as  noise-induced 
transitions  between  different  coherent  patterns  [2],  [1].  The 
study  of  noisy  swarm  dynamics  has  benefited  from  tools  from 
statistical  physics  applied  to  both  first  and  second  order  phase 
transitions  that  have  been  found  in  the  formation  of  coherent 
states  [21]. 

One  important  aspect  of  the  understanding  and  design  of 
space-time  behavior  in  communicating  robotic  systems  is 
that  of  time  delay.  Time  delay  arises  in  latent  communication 
between  agents,  as  well  as  actuation  lag  times  due  to  inertia. 
Time  delays  can  have  interesting  and  surprising  dynamical 
consequences  in  a  system,  such  as  large-scale  synchro¬ 
nization  [22],  [23],  [24],  and  have  been  used  successfully 
for  control  purposes  [25],  [26].  Many  of  the  initial  time- 
delay  studies  focused  on  the  case  of  one  or  a  few  discrete 
time  delays.  Recently,  more  complex  situations  have  been 
considered  such  as  the  case  of  having  several  [27]  and 
random  time-delays  [28],  [29].  Another  interesting  case  is 
that  of  distributed  time  delays,  i.e.  when  the  dynamics  of  the 
system  depends  on  a  continuous  interval  in  its  past  instead 
of  on  a  discrete  instant  [30], 

In  the  case  of  swarming  systems  in  stochastic  environ¬ 
ments,  it  has  been  observed  that  the  introduction  of  a  discrete 
communication  time  delay  induces  a  transition  from  one 
spatio-temporal  pattern  to  another  as  the  time  delay  passes  a 
certain  threshold  [1].  It  was  shown  in  [1]  how  the  complex 
interplay  exists  between  the  attractive  coupling  and  the  time 
delay  in  the  transitions  between  different  spatio-temporal 
patterns  [31],  [32].  Time  delays  in  robotic  systems  have  been 
also  studied  in  the  contexts  of  consensus  estimation  [17]  and 
task  allocation  [18];  in  the  latter,  the  time  delays  originate 
from  the  period  of  time  required  to  switch  between  different 
tasks. 

In  this  paper,  we  consider  a  swarming  model  with  discrete, 
randomly  distributed  time  delays.  We  explicitly  show  how  a 
distribution  of  delay  times  perturb  the  dynamics  from  the 
single  discrete  case  delay  case  analytically.  We  illustrate  the 
dynamical  effects  of  delay  distributions  with  varying  width 
and  show  that  the  system  is  bistable,  and  very  sensitive  to 
choice  of  initial  starting  conditions. 

II.  Swarm  Model 

We  investigate  the  dynamics  of  a  two-dimensional  system 
of  N  identical  self-propelling  agents  that  are  attracted  to  each 
other  in  a  symmetric  manner.  We  consider  the  attraction 
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between  agents  to  occur  in  a  time-delayed  fashion,  due  to 
the  finite  communication  speeds  and  information-processing 
times.  Specifically,  we  focus  on  the  situation  in  which  the 
time-delay  is  nonuniform  across  agents:  there  is  one  time 
delay  for  every  pair  of  agents  %(=  T/i).  for  particles  i 
and  j.  The  time  delays  T,/s  are  time-independent  and  are 
drawn  independently  from  a  random  distribution  p t(t).  The 
swarm  dynamics  are  described  by  the  following  governing 
equations: 

r  i=Vj,  (la) 

v;  =  (l  -  |v,|2)  m  -  ^  E  (r'(0  - TJ(*  ~  %))-  (lb) 
v  y=i 
¥i 


Summing  Eq.  (5)  over  i  and  using  Eq.  (4),  we  get 

R=^l-|R|2-^f|5r;|2jR 

L  (2R-  5r,-  +  |<Sr,-|2)  5r, 

i=  1 

-  R(0  +  4  E  E  (R(f  -  ty)  +  Srj(t  -  T y))  • 

v  i=lj=l 

¥J 

(6) 

We  now  make  some  approximations  on  the  terms  with  the 
double  sums.  For  the  displacements  from  the  center  of  mass, 
we  have 


for  i=l,2...,N.  The  position  and  velocity  of  the  zth  agent 
at  time  t  are  denoted  by  r,  and  V;,  respectively.  Each 
agent  has  self-propulsion  and  frictional  drag  forces  given  by 
the  expression  term  (l  —  |v,|2)  v,-.  The  coupling  constant  a 
measures  the  strength  of  the  attraction  between  agents  and 
the  communication  time  delay  between  particles  i  and  j  is 
given  by  T,y.  Note  that  in  the  absence  of  coupling  agents  tend 
to  move  in  a  straight  line  with  unit  speed  as  time  tends  to 
infinity. 


III.  Mean  Field  Approximation 

We  carry  out  a  mean  field  approximation  of  the  swarming 
system  by  switching  to  particle  coordinates  relative  to  the 
center  of  mass  and  disregarding  the  noise  terms.  The  center 
of  mass  of  the  swarming  system  is  given  by 
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8rj{t  —  i)pT{x)di  =  0, 


8Vj[t  -  r u) 


(7) 


since  Y!J=i  ^r/(l  —  t)  =  0  by  Eq.  (4).  In  passing  from 
the  discrete  to  the  continuous  averaging  above,  we  argue 
as  follows.  The  expression  7^77  Y4L \  i-t;  8 r  / (t  —  Xij )  is  the 
average  of  5r;-(f)  at  the  N  —\  times  t  —  T(  / .  Since  JV  >  1 
and  the  times  T,-y  are  distributed  with  density  p t(t),  this  is 
approximately  equal  to  /0°°  8rj(t  —  x)pT(x)dx. 

Similarly, 
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———^7  ~  [  R(i-t)pt(t)z/t.  (8) 
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Rw  =  ^Er<(0-  (2) 

ly  i—  1 

We  can  decompose  the  position  of  each  particle  into 

r,  =  R  +  5r;,  i=l,  2...  ,N,  (3) 

where  we’ll  have 

f>i(0=0.  (4) 

i=i 

Inserting  Eq.  (3)  into  the  second  order  system  equivalent  to 
Eq.  (1)  and  simplifying  we  get 

R  +  8?i  =  (l  -  |R|2  -  2R  •  Stj  -  |5f,-|2)  (R  +  8rt) 


¥i 


In  a  purely  heuristic  manner,  we  neglect  all  fluctuation 
terms  5r;-(f)  in  the  dynamics  of  the  center  of  mass  and  obtain 
the  mean  field  approximation: 

R  =  (1  -  |R|2)  R  -  a  (  R(f )  -  jf  R(f  -  .  (9) 

where  we  approximated  -^jr-  w  1,  since  we  are  considering 
large  numbers  of  agents. 

IV.  Bifurcations  in  the  Mean  Field  Equation 

The  behavior  of  the  system  in  the  mean  field  approxima¬ 
tion  in  different  regions  of  parameter  space  may  be  better 
understood  by  using  bifurcation  analysis.  This  mathematical 
technique  will  allow  us  to  show  how  the  parameter  plane  of 
coupling  constant  a  and  mean  time  delay  pT  is  divided  into 
regions  with  different  dynamical  behaviors. 

First  we  show  that  Eq.  (9)  has  a  uniformly  translating 
solution  R(f )  =  Ro  +  Vo  •  f ,  where  Ro  and  Vo  are  constant, 
two-dimensional  vectors.  Inserting  the  uniformly  translating 
state  into  Eq.  (9),  we  get 

7*00 

0  =  (l-|Vo|2)V0-fl^  rpt(T)dT  V0, 


(10) 


since  /J°  pz(x)dx  =  1.  Hence,  the  speed  |  Vq  |  of  the  uniformly 
translating  state  must  satisfy 

poo 

|V0|2  =  1  -a  xpz(x)dx=l-apz,  (11) 
Jo 

where  pz  is  the  mean  of  the  pT  distribution.  We  note  that 
the  direction  of  motion  and  starting  point  Ro  are  arbitrary. 

The  other  state  of  interest  is  the  stationary  state  R(f)  =Ro, 
for  an  arbitrary  constant  vector  Ro.  In  the  two-parameter 
space  (a,jttT),  the  hyperbola  apr  =  1  is  in  fact  a  pitchfork 
bifurcation  line  on  which  the  uniformly  translating  states  are 
born  from  the  stationary  state. 

The  linear  stability  of  the  stationary  state  is  determined  by 
the  solutions  to  the  characteristic  equation  of  Eq.  (9): 

@{X)=a(l-  pT(x)e~x'cdx'j-X  +  X2  =  0,  (12) 

and  so  involves  the  Laplace  transform  of  the  distribution  pT. 

In  our  numerical  simulations  of  system  (1),  we  considered 
a  truncated  Gaussian  distribution: 

{b-%)2 

Jfe  2I?  ifr>0  (13) 

0  if  t  <  0, 


where  JV  is  the  normalization  constant.  Note  that  because 
of  the  truncation.  To  and  X\  are  only  approximately  equal 
to  the  mean  and  standard  deviation  of  pz  and  JV  is  only 

approximately  1  /  -y/ 2  7T  Tp . 

We  approximate  the  Laplace  transform  of  the  truncated 
Gaussian  distribution  by  extending  the  integration  range  to 
the  whole  real  line  and  taking  jY  «  yj 2nx2.  In  addition, 
we  approximate  the  mean  and  standard  deviation  of  pT  as 
pz  «  To  and  <7Z  ~  T| ,  respectively.  The  result  is 

[  px(x)e~UdXK,eX^+Xlalil2.  (14) 

Jo 

We  use  the  above  approximation  to  the  Laplace  transform 
of  pz  to  search  for  Hopf  bifurcation  curves  in  the  (a,  pz) 
plane,  by  taking  X  =  ico  in  the  characteristic  equation  (12). 
The  equation  2>{i(o)  =  0  is  equivalent  to: 

a -a2 -uo  =  aeia^-a2(J?'/2,  (15) 


from  which  we  obtain  the  Hopf  curves  parameterized  by  ft) : 


aH((0) 
Pt  h{co) 


(16a) 

n  =  0,1,.... 
(16b) 


In  the  above  expression  for  pT//(ft)),  the  branch  of  tan  in 
(0,  tt)  should  be  used,  since  the  complex  number  on  the 
left  hand  side  of  Eq.  (15)  is  always  on  the  top  half  plane. 
This  family  of  Hopf  curves  labeled  by  n,  together  with  the 
pitchfork  bifurcation  curve  apz  =  1  are  shown  in  Ligure  1, 
for  various  values  of  az. 


Fig.  1.  Hopf  (blue)  and  pitchfork  (red)  branches  in  a  and  jtiT  space. 
The  standard  deviations  of  the  time-delay  distribution  pT  for  the  panels 
(a)  through  (d)  are  0,  0.2,  0.4  and  0.6,  respectively.  Note  the  change  of 
scale  in  the  abscissae. 


When  <JT  =  0,  the  system  exhibits  a  degenerate  point  at 
a  =  1/2,  pz  =  2  (Lig.  1(a)),  where  the  Hopf  bifurcation 
frequency  becomes  zero.  This  is  similar,  but  not  equivalent 
to  a  Bogdanov-Takens  bifurcation  as  is  known  from  previous 
work  [31],  [32].  Since  the  point  on  the  Hopf  curve  in  a  two- 
parameter  bifurcation  plane  occurs  when  the  Hopf  frequency 
becomes  zero,  we  define  this  point  as  a  Zero  Lrequency  Hopf 
(ZLH)  point. 

Lor  az  >  0,  this  ZLH  point  shifts  and  a  second  ZLH  point 
appears  at  a  — >  °°  and  pz  — >  0.  The  location  of  the  two  ZLH 
points  in  the  (a,  pz)  plane  is  given  by  1  /azra)’  where 

aZFH  =  1)1  V7!  ~  °r)- 

When  az  =  0,  the  behavior  of  the  mean  field  in  the  vicinity 
of  the  ZLH  point  is  relatively  well  understood  [32],  [31] 
and  is  as  follows  (see  Lig.  1(a)).  In  the  region  between 
the  pitchfork  and  the  first  member  of  the  Hopf  family,  the 
stationary  state  is  stable.  A  simulation  of  the  full  system  (1) 
with  parameters  in  this  area  reveals  that  indeed  the  center 
of  mass  of  the  agents  comes  to  rest  as  time  progresses 
and  the  particles  spread  themselves  along  a  ring  with  radius 
1  /yfa.  Roughly  half  of  the  particles  move  clockwise  and  the 
other  half  counterclockwise.  Along  the  first  Hopf  curve,  a 
stable  limit  cycle  is  born  and  the  center  of  mass  begins  to 
oscillate  periodically  on  a  circular  orbit.  Below  the  pitchfork 
bifurcation  curve  apz  =  1,  the  translating  state  is  stable. 
Linally,  we  mention  that  there  is  a  region  of  bistability  in  the 
parameter  region  above  the  ZLH  point  (1/2,  2)  between  the 
pitchfork  curve  apz  =  1  and  the  curve  apz  =  2  (not  shown), 
where  the  center  of  mass  can  either  translate  or  rotate.  On 
the  curve  a  pi  =  2  there  is  a  global  bifurcation  where  the 
radius  of  the  orbit  diverges  and  the  limit  cycle  disappears. 

The  above  discussion  helps  us  understand  the  bifurcation 
planes  in  Ligs.  1(b)  through  1(d).  Most  significantly,  we  see 
that  the  parameter  region  where  the  stationary  state  is  stable 
decreases  in  size  as  the  width  of  the  time  delay  distribution 


Fig.  2.  Two  stable  attractors  for  the  swarm  dynamics.  Here  a  —  2,  fir  =  2, 
and  ox  =  0.15.  The  number  of  particles  is  set  to  be  N  =  150.  The  final  state 
shown  for  both  simulations  is  t  =  300.  Panel  (a)  depicts  the  rotating  state 
at  three  snapshots  at  times  t  =297.6,  298.8,  300,  in  red,  green  and  blue, 
respectively.  Panel  (b)  depicts  the  ring  state. 

widens.  Hence  the  system  has  a  higher  tendency  to  behave 
in  an  oscillatory  manner  for  wider  time  delay  distributions. 
This  effect  has  been  corroborated  in  numerical  simulations 
(results  not  shown). 

V.  Numerical  Simulations 

We  analyze  the  dynamics  of  system  (1)  by  solving  the  sys¬ 
tem  of  DDEs  numerically.  We  use  Heun’s  method  together 
with  quadratic  Lagrange  interpolation  to  evaluate  the  time- 
delayed  terms  of  Eqs.  (1).  Overall,  the  numerical  method 
is  second  order  with  respect  to  the  step-size  At.  For  all 
simulations  we  take  the  agents  to  be  uniformly  distributed 
in  a  random  fashion  within  the  unit  box  0  <  x  <  1  and 
0<y  <  1,  and  each  particle  is  initially  at  rest  vj  =  0. 
Moreover,  since  we  are  interested  in  investigating  the  time- 
asymptotic  behavior,  for  all  numerical  experiments  the  time 
of  integration  is  long  enough  to  allow  transients  to  decay. 

In  [31],  [32]  it  was  shown  that  for  the  parameter  set  a  =  2, 
1  =  2  (fixed  delay)  that  the  system  exhibited  a  bistable  set  of 
solutions.  In  the  rotating  state  solution,  all  particles  collapse 
to  a  point  and  that  cluster  of  particles  rotates  around  a  fixed 
center  in  a  circular  orbit.  The  other  possible  stable  solution 
is  a  ring  state,  in  which  all  particles  distribute  themselves 
uniformly  along  a  circle  and  orbit  around  its  center  at  unit 
speed.  Interestingly,  not  all  particles  traverse  the  ring  in 
the  same  direction;  roughly  half  move  clockwise  and  half 
move  anti-clockwise.  We  will  now  examine  these  two  states, 
but  with  random  delays  given  by  the  truncated  Gaussian 
distribution  in  Eq.  (13). 

Figure  2  shows  the  two  final  particle  distributions  after 
transients  in  a  simulation  with  an  initial  state  of  N  =  150 
randomly  placed  particles,  and  where  jXT  =  2  and  ov  = 
0.15.  In  this  case,  depending  upon  the  random  selection  of 
delays,  either  stable  solution  (ring  or  rotating)  is  possible.  To 
understand  the  effects  of  increasing  the  standard  deviation  of 
the  random  delays,  we  use  a  Monte  Carlo  method.  At  100 
different  values  of  (7-  in  the  range  0  <  <7Z  <  0.5,  we  generate 
random  time  delays  from  the  distribution  in  Eq.  (13),  we  then 
simulate  the  system  starting  from  the  same  initial  condition 
and  we  determine  what  state  is  acquired  by  the  swarm  in 
the  long-time  limit.  To  determine  this,  we  first  measure  the 
time-averaged  distance  of  particle  j  to  the  center  of  mass 


Fig.  3.  As  ffT  is  increased,  we  see  a  bifurcation  from  the  stable  rotating 
state.  Panel  (a)  captures  the  transition  from  the  rotating  state  to  the  ring 
state  as  the  standard  deviation  of  the  random  delay  increases.  Panel  (b) 
shows  the  probability  of  converging  to  the  ring  state  for  a  given  aT  of  the 
delays.  These  results  were  compiled  using  a  Monte  Carlo  simulation  with 
100  random  distributions  of  delays  for  100  uniformly-spaced  values  of  oT 
and  for  N  —  50  particles.  See  accompanying  online  movie  and  Appendix  to 
see  the  agents  converge  to  each  stable  pattern. 

over  the  interval  {t\ .  ti)'. 

(8rj)(tut2)  =  r~T  i  l<5r/(f)l<*  (!7) 

where  the  size  of  the  interval  (t i,  12)  is  long  enough  to 
include  several  periods  of  oscillation.  The  ensemble  average 
of  Eq.  (17)  is  then 

(Sr)(h,t2)  =  i}L(5ri)(tut2y  (18> 

j= 1 

A  value  (5r)(?]  ,2)  ~  1  /^/a  will  indicate1  that  the  system  has 
converged  to  the  ring  state,  while  (tSr)/^  tl\  -C  1 1 \fct  shows 
that  the  rotating  state  has  been  adopted  instead2. 

Figure  3(a)  demonstrates  the  effect  of  increasing  az  on 
the  final  state.  The  blue  circles  show  that  for  <7r  small,  this 
initial  condition  converges  to  the  rotating  state.  However,  for 
<7Z  >  0.2  the  same  initial  collection  of  particles  will  converge 
to  the  ring  state  with  high  probability.  In  between,  there  is  a 
transition  region  where  both  states  are  commonly  observed; 
the  state  that  occurs  depends  on  the  random  choice  of  time 
delays.  The  black  dashed  lines  of  3(a)  show  two  simulations, 
one  which  starts  near  the  rotating  state  (the  lower  curve),  and 
one  which  starts  near  the  ring  state  (the  upper  curve)  as  <7T 
is  increased.  These  curves  demonstrate  the  stability  of  these 
steady  states,  and  the  effect  of  random  delays  near  these 
states. 

Figure  3(b)  shows  the  conditional  probability  of  ending 
up  in  the  ring  state  as  a  function  of  (7-.  As  expected,  for 
this  choice  of  initial  conditions,  for  <7T  small  enough,  there 
is  zero  probability  of  leaving  the  rotating  state;  however,  as 
crT  is  increased,  the  probability  increases  to  one. 

The  results  of  these  numerical  studies  strongly  suggest 
that  even  though  there  is  bi-stability  between  the  ring  and 
the  rotating  states,  the  size  of  their  respective  basins  of 
attraction  is  changing  dramatically  as  the  standard  deviation 
cjt  increases. 

'When  the  delays  are  uniform,  the  ring  state  has  a  radius  of  1  j\fa  [32], 

2This  is  true  for  the  range  of  values  of  aT  considered  here. 


VI.  Discussion 

In  this  paper  we  studied  the  dynamics  of  a  self-propelling 
swarm  with  time-delayed  inter-agent  attraction.  In  contrast  to 
the  previously  considered  case  of  uniform  time  delay  across 
agents,  we  considered  the  situation  in  which  the  time  delay 
between  every  pair  of  agents  is  drawn  randomly  from  a 
distribution  pT. 

Using  a  mean-field  model  of  the  swarm,  we  showed  how 
the  two  parameter  bifurcation  plane  of  coupling  strength 
and  mean  time  delay  changes  with  respect  to  the  case  in 
which  all  time  delays  are  equal.  The  full  implications  of 
these  bifurcation  results  are  the  subject  of  our  ongoing  work. 
In  particular,  it  is  unclear  what  the  stable  solutions  are. 
Nevertheless,  the  dramatic  changes  seen  in  the  two  parameter 
bifurcation  plane  as  the  standard  deviation  crr  increases 
suggest  that  the  basins  of  attraction  of  each  attractor  undergo 
big  changes  as  well. 

Our  numerical  experiments  show  that  the  swarm  displays 
bi-stable  behavior  between  the  ring  and  rotating  states,  at 
the  parameters  considered.  Interestingly,  however,  our  work 
suggests  that  the  basin  of  attraction  of  the  ring  state  greatly 
expands  as  the  distribution  of  time-delays  pT  widens.  Thus, 
in  a  sense,  widening  the  distribution  of  time-delays  stabilizes 
the  stationary  state  of  the  swarm  center  of  mass. 

Even  though  in  our  model  the  attractive  force  among 
agents  is  linear,  we  believe  this  work  is  useful  since  it 
represents  a  first  approximation  for  other,  more  general  forms 
of  attractive  interaction.  Here,  we  have  limited  our  focus  to 
the  case  where  the  delays  between  agents  are  symmetric  and 
constant.  However,  one  important  generalization  of  this  sys¬ 
tem  involves  incorporating  time  dependent  delays,  including 
those  which  vary  as  a  function  of  the  distance  between  the 
two  agents.  This  particular  refinement  of  our  model  is  the 
subject  of  ongoing  work  and  beyond  the  scope  of  the  current 
paper. 

Finally,  although  we  did  not  consider  repulsion  between 
agents,  preliminary  research  leads  us  to  believe  that  the 
patterns  observed  in  this  investigation  persist  when  the 


characteristic  repulsion  strength  between  robots  is  small 
compared  to  global  attraction  parameters.  For  these  reasons, 
our  results  indicate  how  to  exploit  time-delayed  actuation 
when  designing  swarm  robotic  systems  with  desired  tasks 
and  functionalities. 

VII.  Appendix-Video  Description 

The  purpose  of  this  research  is  to  investigate  the  effects 
of  randomized  communication  delay  on  emerging  patterns  in 
swarming  dynamics.  This  short  video  captures  the  transition 
between  two  different  stable  patterns  for  a  swarm  as  a 
function  of  the  standard  distribution  of  the  delays. 

The  two  coordinate  axes  in  the  video  show  a  scatter  plot 
of  the  positions  of  the  particles  animated  in  time.  The  initial 
positions  are  identically  randomly  distributed  particles  in  the 
unit  box.  The  temporal  state  of  the  swarm  is  updated  in 
time  using  a  numerical  scheme  called  Heun’s  Method,  and  a 
snapshot  is  captured  at  every  discrete  time  interval.  Here,  the 
left  coordinate  axis  uses  a  standard  deviation  of  the  delays 
crx  =0.1  while  the  right  axis  uses  a  standard  deviation  of 
crT  =  0.3.  The  mean  delay  for  each  simulation  is  set  to  be 
jj.z  =  2,  and  the  number  of  particles  for  both  is  N  =  50. 
The  vectors  at  each  particle  give  the  velocity  associated 
with  that  particle.  The  two  simulations  are  run  side  by  side 
to  demonstrate  the  dynamics  involved  in  converging  to  the 
“rotating”  final  state  on  the  left,  and  the  “ring”  final  state 
on  the  right.  The  video  demonstrates  the  dynamics  over  the 
time  interval  from  t  =  0  to  t  =  45,  and  so  includes  transients. 
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